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Abstract 

We present the crossover line between the quark gluon plasma and the hadron 
gas phases for small real chemical potentials. First we determine the effect of 
imaginary values of the chemical potential on the transition temperature using 
lattice QCD simulations. Then we use various formulas to perform an analytic 
continuation to real values of the baryo-chemical potential. Our data set main¬ 
tains strangeness neutrality to match the conditions of heavy ion physics. The 
systematic errors are under control up to /r_B « 300 MeV. For the curvature of 
the transition line we find that there is an approximate agreement between val¬ 
ues from three different observables: the chiral susceptibility, chiral condensate 
and strange quark susceptibility. The continuum extrapolation is based on Nt = 
10, 12 and 16 lattices. By combining the analysis for these three observables we 
find, for the curvature, the value k = 0.0149 ± 0.0021. 


1. Introduction 

For heavy ion physics, the most important feature of the phase diagram 
of Quantum Chromodynamics (QCD) is the line that separates the hadron gas 
phase from the quark gluon plasma, and the conjectured critical end-point along 
this line separating cross-over from first order transition [T]. 

The qualitative form of the phase diagram was sketched four decades ago 
[2] as a consequence of Hagedorn’s exponential spectrum of hadron masses [3] . 
The order of the transition at zero density has been determined much later, 
for Nature’s selection of quark masses the two high temperature phases are 
connected through a cross-over [4] . In the absence of a real transition the cross¬ 
over temperature can be determined but it is ambiguous [5]. Observables that 
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are related to the spontaneous breaking of chiral symmetry (chiral condensate 
and its susceptibility) give a temperature around 155 MeV [SHE]. 

Beyond the transition temperature Tc at vanishing density, the chiral cross¬ 
over line is described by a standard curvature parameter (k) and higher order 
terms: 



( 1 ) 


Extracting Tc{^,b) from first principles is very challenging. Direct Monte- 
Carlo calculations are hindered by the sign problem. Attempts to reach non¬ 
vanishing fj,B on the lattice include reweighting of the generated configurations 
[HHS], Taylor expansion in fi [HHiH], analytic continuation from imaginary fj, 
, use of the canonical ensemble and density of state methods [sni 

[8T] . More recent approaches are represented by the use of dual variables [32] , and 
the complex Langevin equation [831134] . However, their application to QCD with 
physical parameters and controlled discretization has not yet been achieved. 
The phase diagram was frequently studied in various model frameworks, see 
e.g. Ref. [T] and references therein. Recently, functional methods have also 
been applied to QCD |3SH37]. 

For the first few coefficients in Eq. Q it is enough to study QCD at small 
fiB-: for which there are several methods, k can be and has been determined 
by calculating the ^s-derivative of the chiral condensate using only ^b = 0 
ensembles [IH1I3HI- However, the signal/noise ratio of higher ^b derivatives is 
suppressed with powers of the volume, making this approach impractical beyond 
order. Lattice calculations are perfectly feasible, though, with imaginary 
values of the chemical potential [T0l|T9l|39]. Setting ^b = ifJ-s avoids the 
sign problem and the transition line can be studied |iDM5] . 

In this study we follow the imaginary-p.^ approach and go beyond previ¬ 
ous studies by a) performing a continuum approximation with lattices up to 
Nt = 16; b) tuning pls{^b,T) such that the strangeness neutrality condition 
is maintained; c) using several observables: chiral condensate, chiral suscepti¬ 
bility and strange susceptibility; d) comparing the Taylor and the imaginary-/^ 
method for the strange susceptibility; e) calculating the systematic errors from 
scale setting, fit ranges, analytic formulas, etcQ 

This Letter first gives a brief account of the necessary lattice simulations at 
zero and finite temperatures. Then the method for setting strangeness neutrality 
is explained. Finally, we give a detailed description of the analysis and present 
the continuum results for the curvature. 


^During the writing of this manuscript a similar independent analysis, based on Nt = 
6, 8,10,12 lattices and analytic continuation from imaginary fig appeared in arXiv 1441 . Their 
findings are similar to ours but the present analysis has finer lattices, smaller pion splittings 
and significantly larger statistics. 
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2. Simulation setup 

This study is part of the 2nd generation staggered thermodynamics program 
of the Wuppertal-Budapest collaboration [iS]. We use a four times stout [15] 
smeared {p = 0.125) staggered fermion action with 2+1+1 flavors, i.e. dynam¬ 
ical up, down, strange and charm quarks. The gauge action uses the tree-level 
Symanzik improvement. The two light quarks are degenerate, their masses and 
the strange quark mass are tuned such that the physical pion and kaon mass 
over pion decay constant are reproduced for every lattice spacing. For the zero 
temperature runs we kept the volume large LtUtt > 4 in the entire lattice spac¬ 
ing range of interest for this study: a = 0.2 ... 0.063 fm. The charm mass was 
set to rric/rns = 11.85 [47]. The simulation parameters are detailed in Ref [45] . 
The overall scale was determined from /^. We used Wq as an alternative scale 
setting for the analysis [48] . 

The chiral susceptibility as well as the chiral condensate require renormal¬ 
ization. The additive divergence is removed by subtracting the vacuum expec¬ 
tation value, the multiplicative divergence canceled by the same factor in the 
bare quark mass mi- The renormalized condensate and its susceptibility are 
dimensionful quantities, we use the fourth power of the pion mass to form a 
dimensionless observable. We do not restrict the chiral susceptibility to the dis¬ 
connected part. The third observable that we use to identify the /i-dependent 
transition temperature is the strange susceptibility: thanks to the exact quark 
number conservation it does not require renormalization. 

At finite temperature, we have collected data at zero and at imaginary baryo- 
chemical potentials. The ps = 0 data are used to perform a Taylor expansion 
on one of our studied observables, and also to obtain a “baseline” for the shifted 
transition temperatures at other ps values. The zero density configurations are 
listed in Ref. [45] . 

The range of imaginary baryo-chemical potentials is limited by the Roberge- 
Weiss transition at ps = *7rT [33]. Below a limiting temperature there is 
no transition as Im [pb/T] crosses tt, but there is a first order transition above 
Trw, where the imaginary density is non-vanishing and flips sign at Pb/'^ ~ 
The nature of the transition at T^w depends on the quark masses [501152] . For 
intermediate masses the system at T = Tj^w ^nd Pb/'^ ~ critical, 

and then in the entire range of smaller imaginary chemical potentials we will 
see a crossover in temperature. Our data suggests that, for physical masses, the 
latter scenario is realized, namely we are working with a cross-over for all used 
Ms/r. 

We selected six imaginary chemical potential values: 

i = 1,2,3,4, 5, 6 (2) 

We have all six j values for our Nt =8, 10 and 12 lattices and only j = 3... 6 
for Nt =16. The reason for this is the following: j = 0... 5 data are needed to 
determine the simulation parameters at finite imaginary such that the 
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strangeness neutrality condition is fulfilled (see later). The continuum extrap¬ 
olation for this analysis can be carried out using Nt =8, 10 and 12 lattices. 
For the determination of the curvature k of the phase diagram we also need 
the finest Nt = 16 lattices. Since j =1 and 2 do not give a statistically very 
significant contribution to k we decided not to have these two points in our most 
expensive Nt = 16 ensembles. Therefore, in order to have the same setup for 
all lattice spacings, the k determination is based on j = 3,4 and 5. The j = 6 
point is used to estimate higher order effects. 

This range to find the k coefficient (l^g/T < 2) is narrower than in ear¬ 
lier studies (e.g. < 2.36 in Ref. [H] and < 2.6 in Ref. [45]1. 

A broader range of chemical potentials has the advantage that the numerical 
derivative \Tc{^b) — 7c(0)]/Mb ^ larger signal/noise ratio. However, more 

non-linearities appear in a broader range and the results are more prone to sys¬ 
tematic errors as the singularity at « ttT is approached. This is the reason 
(to avoid unwanted systematic uncertainties) why we have taken a smaller 
range and we use other methods to increase the signal/noise ratio. 

We performed simulations on 32^ x 8, 40^ x 10, 48^ x 12 and 64^ x 16 
lattices, at sixteen temperatures in the temperature range 135... 210 MeV. We 
have generated between 10000-15000 Hybrid Monte Carlo updates, analyzing 
every 5th of them (every 10th for Nt = 16). The configurations have been 
evaluated for up to fourth order generalized quark number susceptibilities [53] 
and for the chiral condensate and susceptibility. For = 0 we have 5... 10 
times more statistics, this ensures a solid guidance to the fitting procedure. 

3. Strangeness neutrality 

The most popular representation of the QCD phase diagram is in the tem¬ 
perature vs. chemical potential plane. The baryo-chemical potential axis leaves 
room for various interpretations. Ref. [42] used the full baryo-chemical poten¬ 
tial including the strange quarks, i.e. = fJ-d = fJ-s = Mb/ 3- In Ref. |33| both 
= (J^b/^ and = Md = Mb/3, Ms = 0 were studied. 

However, neither of the recipes Ms = 0 or = /iB/3 maps consistently to 
the situation that is realized in experiment. In heavy ion collisions, non-strange 
particles are colliding. Although ss are generated in the collision, the net- 
strangeness is zero. Therefore we want to tune the chemical potentials to such 
values which guarantee strangeness neutrality. The light chemical potentials 
are kept identical (mu = Md), which ensures isospin symmetry also at finite mb- 
This corresponds to an experimental situation where Z = 0.5A. Alternatively, 
one can achieve a different ZjA ratio corresponding to heavy nuclei by tuning 
Mu and Md appropriately. This possibility will be discussed later. Requiring 
strangeness neutrality and fixing the value of ZjA (or alternatively the electric 
charge/baryon number ratio) uniquely determines all three quark chemical po¬ 
tentials as functions of mb- For the isospin symmetric case m« = Md = Mb/3 
so the only non-trivial task is to find the strange quark or strangeness chemical 
potential. 
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The strange quark chemical potential (fig) is related to the strangeness (^ 5 ) 
and baryo-chemical potential (/is) as fig = — /is- Then /ig = /xs/S 

approximates strangeness neutrality at low temperature, and /Xg = 0 at high 
temperature. In this work we determine the strangeness neutral trajectory 
from lattice simulations. 

It is relatively straightforward to perform Taylor expansions from /x^ = 0 
on the trajectory that respects strangeness neutrality. For the equation of state 
|54[ 155) and for fluctuations relevant for calculating freeze-out parameters in 
heavy ion collisions [SH [56] this procedure is already standard. 

For actual simulations at finite /x^ the strange chemical potential has to be 
fine tuned for every temperature, baryo-chemical potential and lattice spacing. 
We solved this challenge by solving the 


d dlogZ 
d/x^ dus 


= 0 , 


(3) 


differential equation discretized in /x^ with the trivial initial condition d log Z/dfis 
0 at = 0. This equation simply states that the derivative of strangeness 
is zero. Using the 2nd order explicit Runge-Kutta scheme, we determine iJ-sids) 
using the prescription: 


dsid-B + — d-sid-B ~ ~ 2 

Xs 


Afi 


B I 


(4) 


with the step size A/x^/T = tt/S (see Eq. [^. For the initial step (/x|j/r = 
Afig/T) we used the high-statistics /xg = 0 runs and NLO Taylor expan¬ 
sion. Each step using Eq. Q requires a simulation at /x^ and the evalua¬ 
tion of the 2 nd order fluctuations: Xs^b — ^/{TV)d‘^ log Z/dfisdfJ,B and yf = 
l/{TV)d'^ log Z/9/x|. This method would be 0(A/x^ ) accurate only, but as an 
additional correction, we do a small extrapolation for both terms on the RHS 
of Eq. Q after each simulation so that the remaining strangeness neutrality 
violation is not propagated to the next step. For this extrapolation, we need 
higher order fluctuations [53]. This combination is 0(A/x^^) accurate in the 
complete /x^ range. The resulting fis{dB,T) function is interpolated in T and 
extrapolated in l/N^ and the resulting smooth function is used to start the 
simulations at -I- A^^. In Fig. [^we show the resulting strangeness chemical 
potential. 


4. Analysis details 

We calculate the curvature of the phase diagram from three observables. We 
calculate statistical and systematic errors for all three. 

1) Our first observable is the chiral susceptibility Xipip/T^t- discussed 
previously, it requires additive and multiplicative renormalization. For details 
on this procedure see Ref. |7]. The chiral susceptibility forms a peak at the 
transition temperature. With increasing imaginary chemical potential this peak 
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Figure 1; The imaginary strangeness chemical potential that realizes strangeness neutrality. 
Here we show a continuum extrapolation based on 32® X 8, 40® X 10, 48® X 12 lattices for 
/Tg/r = 0.39 and 0.79 and also using 64® X 16 for the larger chemical potentials. The error 
is statistical only. 


is shifted towards higher temperatures, approximately maintaining its height 
and width. For other normalizations (e.g. the shape of the function 

changes more significantly while varying the chemical potential. 

We fit T)/m^ in a global fit function where for each /is a different 

width, height and peak position is allowed, but the other parameters that de¬ 
scribe the peak shape are /is-independent. We use two different modifications 
to the Lorentzian peak form: 

^ic + A^p){l + W\tJ.){T-T,{^i))^y forT<r, 
mt \c + A^{p) ll + B^W\^i){T-T,iii))y for T > T, 

and 

= r* 4- __ _ ('c'l 

The dependent parameters A{fi), W{fi) and Tc{^j,) describe the change in 
the height, width and the position of the curve as ^ increases. For the zero 
temperature data which are required for renormalization we use two different 
interpolations in the inverse gauge coupling: a 6th order polynomial and a 
simple rational function. We have two options for the scale setting using /^r or 
Wo and we apply three possible fit windows to select the transition range. In 
order not to interfere with the shifted temperature dependence the fit windows 
constrain the value of the susceptibility, not the temperature. 

The effect of the /i dependent parameters is a shift in T, and a rescaling in 
T and y. Applying the inverse transformation to the finite /ig data points all of 
them should collapse on the = 0 curve. This is demonstrated in Fig. The 
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Figure 2: The chiral susceptibility at several imaginary chemical potentials on our 48^ x 
12 lattice. After a /i^-dependent shifting and stretching, data from all chemical potentials 
collapse on one curve. The fitted curve corresponds to eqn. § at fiB = 0. 


advantage of this procedure is that the /r independent parameters can be fitted 
using the the high-statistics runs at = 0 and the non-vanishing runs are 
needed to determine the relative position and rescaling compared to this more 
complicated functional form. This allows the precise determination of AT'c(/i^) 
with an error below 0.25 MeV, while Tc{hb) itself has an error of several MeV. 
We extract k by a linear fit of ATc vs. in the range 1.2 < < 2, and 

extrapolate k to the continuum. Since the continuum extrapolation of k had 
a large when all four lattices were used we included only Nt = 10,12 and 
16 in the final result, resulting in a good for analyses. In an alternative 
analysis we made a combined continuum and fit again using only the finest 
three lattice spacing, and found acceptable x^ values again. 

2) The chiral condensate {'^4’Y = rnq{d\og Z/dmq )/is a remnant order 
parameter of the chiral transition. Its inflection point (though it is hard to 
locate in a finite precision data set) is very close to the peak position of 
At finite Hg the temperature dependence of {'tpYY is shifted and very slightly 
stretched. 

We find that the data at fiB = 0 (see Ref. [7]) can be very accurately 
described by simple fit functions. The /x dependence in this case is well described 
by just two /is-dependent parameters describing a shifting and rescaling of the 
renormalized condensate. We use the following parameterizations: 

T) = A{fi) {1 + B tanh [C {T - T^Y)] + D {T - TYfi))) (7) 

and 

T) = A{p.) {l + B arctan [C (T - + D {T - r,(/x))). (8) 

Similarly to the chiral susceptibility, we use two possible zero temperature in¬ 
terpolations (6th and 7th order polynomials of the inverse gauge coupling), two 
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scale settings, four fit windows, k is obtained either from a combined and 
continuum fit, or separately. 

3) The analysis of the strange susceptibility yf go®® along the lines of the 
chiral condensate. A simplification here is the absence of renormalization. Since 
this quantity is the most sensitive to the actual value of strangeness, before the 
analysis we correct for the inaccuracies of the strangeness neutrality condition 
using the higher ns fluctuations. Although its inflection point does not have to 
agree with that of the chiral condensate, we find that the shifting effect of the 
chemical potential is very similar. 

For all three quantities we make a histogram of the results from all analyses. 
For the chiral susceptibility we have two T > 0 fit forms, two T = 0 interpola¬ 
tions, two scale settings, three fit windows and either separate or combined k 
extraction and continuum limit. This results in 2 • 2 • 2 • 3 • 2 = 48 analyses. For 
the condensate we have the same choices but with four fit windows resulting in 
64 analyses. For the strange susceptibility there is no renormalization, thus no 
T = 0 interpolation is needed which leads to 32 analyses. The central 68% of 
the histograms estimates our systematic error. The statistical error is obtained 
from 1000 bootstrap samples. The two errors are of similar magnitude and they 
are added in quadrature resulting in our final uncertainties. 

We summarize our results for the curvature in Table [TJ 


Chiral susceptibility 
Chiral condensate 
Strange susceptibility 
Susceptibility at Z = 0.4A 


0.0158 ±0.0013 
0.0138 ±0.0011 
0.0149 ±0.0021 
0.0149 ±0.0017 


Table 1: The curvature (k) of the QCD phase diagram in the continuum limit from various 
observables, k is fitted in the range 1.2 < << 2. 


The histograms of the three quantities can be joined into a single one lead¬ 
ing to our combined result based on our three observables with strangeness 
neutrality: 

K = 0.0149 ± 0.0021. (9) 

We also consider the curvature for the case when not only the strangeness 
neutrality, but also proper charge/baryon density ratio is reproduced (for lead 
and gold ions: Z « 0.4A). We achieve this by Taylor-extrapolating the strange 
susceptibility for every finite ns ensemble to leading order, and fitting as be¬ 
fore. We conclude that the difference between Z = 0.4A and Z = 0.5A phase 
diagrams is negligible for small ns- 

For small enough imaginary chemical potential the analytical and the Taylor 
method have to give the same curvature at every lattice resolution. In the 
Taylor method one expands the observables in fj-B, the leading coefficients are 
calculated from = 0 simulations and then used to extrapolate to finite hb- 
Fig.j^shows how this expansion compares to the direct simulations for our j = 5 
chemical potential which is the largest one used to extract k. A comparison in 
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Figure 3: Comparison of the strange susceptibility obtained from direct simulation and 

extrapolation for = 57r/8. The blue circles and squares correspond to the strangeness 

neutral case obtained via extrapolation from // = 0 and direct simulations, respectively. The 
green triangles show the full baryo-chemical potential case obtained via extrapolation from /i = 
0. There are no direct simulations in this case but one can extrapolate from the strangeness 
neutral direct point (green dots). As a reference the /.t = 0 data are also shown (red crosses). 


the case of full baryo-chemical potential is also shown. The extrapolated data 
are then fitted for k as if they were simulated at finite At Nt = 10 we find 
K, = 0.0131(9) from the direct simulations and k = 0.0115(10) from the Taylor 
expansion. The agreement indicates that we are still in the linear regime and 
the extraction of k using j = 3,4, 5 is safe. 

Finally we estimate the systematics of the extrapolation to real fj,B- We 
include the j = 6 data points into the analysis and allow for non-linear 
fits. We consider fitting Tc{fi%)/Tc with the functions 1 -I- ax, \ + ax + bx"^, 
(1 -I- aa;)/(l -I- bx) and [1 + ax + bx)~^ with x = All these functions are 

analytic in x and they represent various analytic continuations of the Tc{fig) 
imaginary ^ phase diagram. The difference between these ansatzes provides a 
systematic uncertainty for the real fj, phase diagram. 

Our main results are depicted in Fig. Since the curvature from both the 
strange susceptibility and from the chiral condensate/susceptibility are consis¬ 
tent with each other we show only one curve. The curvature from the chiral 
condensate is our most precise result, therefore we present the transition line 
coming from this observable. The corresponding transition temperature at fi=0 
is at 157 MeV. At intermediate real we observe a significant rise in the un¬ 
certainty due to the statistical error on the non-linear /i^-dependence and the 
ambiguity of the analytic ansatz. This sets the range of validity for this study. 

The present result indicates a stronger curvature than the one presented 
in Ref. [ 35 ] • There are, however a couple differences between the definitions/ 
approaches of the curvature of the present analysis and Ref. [ 35 ]. Note that 
the transition is a smooth cross-over, thus different definitions obviously lead to 
different results. 
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Figure 4: The phase diagram based on the /^-dependent Tc from the chiral condensate, 
analytically continued from imaginary chemical potential. The blue band indicates the width 
of the transition. The shaded black region shows the transition line obtained from the chiral 
condensate. The widening around 300 MeV is coming from the uncertainty of the curvature 
and from the contribution of higher order terms, thus the application range of the results 
is restricted for smaller fi values. For completeness, on the right panel we also show some 
selected non-lattice results: the Dyson-Schwinger result of Ref. m and the freeze-out data 
of Refs. |57H63| . 


a. In Ref. [35] we used a vanishing strangeness chemical potential. In the 
present analysis we use instead vanishing strange density. The reason for this 
change is to be as close to the experimental situation as possible. In heavy ion 
collisions the net strangeness is zero. 

b. It is emphasized in the discussion of Figure 5 of [35] that only statistical 
uncertainties were provided. The present analysis estimates systematic uncer¬ 
tainties coming from various aspects of the analysis as discussed earlier. These 
are comparable to or in some cases even larger than the statistical uncertainties. 
A similar assumption on the systematics of Ref. |38j would make the tension 
between the results much weaker. 
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